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Abstract 

Numerical  efficiency  comparisons  of  a  four-node  finite  element  model  (FEM),  a 
mass-spring  lattice  model  (MSLM),  and  a  mass-spring-dashpot  lattice  model  (MSDLM)  are 
investigated.  Specifically,  the  error  in  the  ultrasonic  phase  speed  with  variations  in  Poisson’s 
ratio  and  angle  of  incidence  is  evaluated  in  each  model  of  an  isotropic  elastic  solid.  With  regard 
to  phase  speed,  materials  with  constant  N  grid  spaces  per  P- wavelength  having  Poisson’s  ratios 
between  0.0  and  0.25  are  modeled  more  accurately  with  the  MSLM.  Materials  with  Poisson’s 
ratios  between  0.35  and  0.5  and  N  grid  spaces  per  P-wavelength  are  more  accurately  modeled 
with  the  FEM.  Materials  whose  Poisson’s  ratio  is  between  0.25  and  0.35  are  modeled  equally 
accurately.  With  regard  to  phase  speed,  viscoelastic  materials  modeled  with  FEM  and  MSDLM 
show  good  agreement  with  known  analytical  solutions.  The  computational  expense  of  all  three 
models  is  also  examined.  The  number  of  floating  point  operations  (FLOPS)  needed  to  achieve  a 
specified  phase  speed  accuracy  is  calculated  for  each  different  model.  While  the  FEM  and 
MSLM  have  nearly  the  same  computation  cost,  the  MSDLM  is  5  times  more  costly  than  either 
the  FEM  or  MSLM. 
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Introduction 

Nondestructive  evaluation  (NDE)  techniques  have  been  used  for  decades  to  characterize 
materials  and  inspect  products.  For  example,  the  U.S.  Navy  has  developed  a  variety  of  NDE 
techniques  and  systems  for  identifying  defects  in  ship  structures.  These  techniques  and  systems 
are  vital  for  finding  defects  before  they  impact  the  safety  and  mission  readiness  of  a  ship  [1]. 
Most  NDE  systems  contain  an  energy  source  used  to  probe  an  object,  a  receiver  or  detector  that 
measures  how  the  energy  has  been  changed  by  the  object,  and  components  and  analyses  to 
record,  process,  and  interpret  the  measurement  data  [2].  Many  of  these  systems  use  ultrasonic 
wave  energy.  Some  advantages  of  ultrasonic  testing  (UT)  are  that  small  surface  and  subsurface 
discontinuities  can  be  detected.  Ideally,  the  approximate  size  and  orientation  of  the  flaw  can  also 
be  detennined  [3].  Prediction  and  simulation  of  ultrasonic  wave  propagation  provide  valuable 
analytical  techniques  in  the  interpretation  of  UT  data. 

Two  approaches  to  numerical  simulation  of  wave  propagation  are  with  a  finite  element 
model  (FEM)  and  a  mass-spring-dashpot  lattice  model  (MSDLM).  Finite  element  procedures  are 
now  a  critical  part  of  engineering  analysis  and  design  [4].  The  versatility  and  ease  of  use  of 
commercially  available  finite  element  programs  have  contributed  to  their  popularity.  In  recent 
years,  new  research  in  mass-spring  lattice  models  has  taken  place.  Yim  [5]  discusses  the 
advantages  of  a  mass-spring  lattice  model  and  the  U.S.  Navy  has  already  begun  research  in 
numerical  simulation  of  thick,  layered  composites  with  mass-spring-dashpot  lattice  models  [6]. 
In  this  thesis,  FEM  and  MSDLM  are  compared  with  regard  to  phase  speed  accuracy  and 
computational  cost. 
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Problem  Statement 

Choosing  an  appropriate  analysis  method  can  be  difficult  without  knowing  which  model 
will  most  accurately  represent  a  particular  material  or  phenomenon. 

Numerical  anisotropy  exists  in  each  model  and  may  affect  computational  results, 
especially  with  regard  to  propagation  direction  [7].  It  is  well  known  that  one  method  of  reducing 
error  is  by  reducing  mesh  size.  This  does  have  the  added  effect  of  increasing  computational  cost. 
This  thesis  proposes  that  by  choosing  carefully  the  model  used,  error  and  computational  cost 
may  be  reduced. 

Introduction  to  Numerical  Models 

The  mass-spring-dashpot  lattice  model  (MSDLM)  in  Fig.  1  [8]  is  a  modified  version  of 
the  mass-spring  lattice  model  (MSLM)  of  Yim  and  Sohn  [9]  and  Yim  and  Choi  [10].  The 
MSDLM  is  a  heuristic,  physical  model  where  the  inertia  and  viscoelasticity  of  a  solid  are 
modeled  as  particles  interconnected  with  springs  and  dashpots. 

The  spring  constants  and  dashpot  coefficients  are  derived  from  the  exact  partial 
differential  equations  governing  a  two-dimensional  standard  linear  solid  [8].  The  particle 
velocities  and  displacements,  as  well  as  volumetric  forces  through  each  element  are  numerically 
integrated  with  a  fourth-order  Runge-Kutta  algorithm  [11].  The  MSDLM  has  been  used  to  study 
wave  propagation  phenomena  in  materials  having  attenuation  and  has  been  shown  to  agree  with 
analytical  solutions  in  both  steady-state  and  transient  analyses  [8], 
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Figure  1.  Schematic  of  an  MSDLM  at  an  interior  plane-strain  particle  located  at  position  (i,j)  [8], 

To  ensure  stability  and  convergence  of  the  Runge-Kutta  algorithm  in  the  MSDLM,  the  Courant 
number  C  must  satisfy 
c  „A  t 

c_  max,P  <  1  3Q  (i) 

h 

where  cmax,p  is  the  maximum  P  phase  speed,  h  is  the  grid  space,  and  At  is  the  numerical  time 
step  [8],  Phase  speed  is  the  speed  at  which  the  crest  of  a  single-frequency  wave  travels. 

The  finite  element  method  is  an  approach  for  solving  partial  differential  equations  (PDEs) 
and  integral  equations  [12].  Finite  element  modeling  of  a  solution  involves  the  mesh 
discretization  of  a  continuous  domain  into  a  set  of  discrete  sub-domains  and  a  finite  number  of 
points  called  nodes.  In  Fig.  2,  elements  comprising  the  entire  domain  are  connected  at  common 
nodes  and  collectively  approximate  the  shape  of  the  initial  domain.  The  elements 
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are  then  approximated  by  a  discrete  set  of  piecewise  continuous  functions  (polynomials)  defined 
using  the  nodal  values  of  the  continuous  solution.  The  solution  to  the  piecewise  continuous 
functions  approximates  the  solution  to  the  initial  PDE.  Finite  difference  methods  are  different 
from  finite  element  methods  in  that  the  differential  operators  are  approximated  rather  than 
represented  by  piecewise  polynomials. 


i-lj-1  i,j+l  i+l,j+l 


Figure  2.  Schematic  of  a  four-node  FEM  at  an  interior  plane-strain  node  located  at  position  (i,j). 

In  structural  mechanics,  finite  element  methods  are  often  based  on  an  energy  principle 
such  as  the  principle  of  virtual  work.  Polynomial  shape  functions  are  used  to  relate  the 
displacement  at  any  particular  point  to  the  displacements  at  the  FEM  nodes.  A  model  of  a 
material  such  as  a  standard  linear  solid  can  be  generated  using  the  displacements  at  the  nodes, 
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material  properties,  and  constitutive  relationships.  FEM  have  been  used  to  study  wave 
propagation  and  attenuation  and  have  been  shown  to  agree  with  analytical  solutions. 

To  ensure  accuracy  and  stability  in  a  wave  propagation  finite  element  model,  certain 
criteria  must  be  met.  The  length,  Le, of  a  finite  element  must  be 

L  =  c  At  (2) 

e  max  V  / 


where  c  is  the  maximum  wave  speed  and  At  is  the  corresponding  timestep.  For  this  analysis,  an 
explicit  integration  method  is  used  which  means  the  Courant-Friedrichs-Lewy  (CFL)  condition 
must  be  met  in  order  to  obtain  accurate  results.  Simply  stated,  the  maximum  allowable  Courant 
number  C  that  an  explicit  time-integrator  may  use  is  1 . 


C  = 


c  At 

max 


(3) 


where  Le  corresponds  to  the  length  of  the  finite  element  h.  While  it  is  well  known  that  C  greater 
than  1  causes  a  model  to  become  unstable,  a  finite  element  model  also  becomes  inaccurate  as  C 
decreases  to  values  less  than  1  [12]. 

Fig. 3  is  a  mass-spring-lattice  model  (MSLM).  It  is  similar  to  the  MSDLM  in  arrangement 
of  nodes  but  differs  from  the  MSDLM  in  that  it  lacks  dashpots.  It  is  useful  in  modeling  wave 
propagation  phenomena  in  elastic  materials  and  has  been  shown  to  agree  with  analytical 
solutions  in  both  steady-state  and  transient  analyses  [5,9,10]. 
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Figure  3.  Schematic  of  MSLM  at  an  interior  plane-strain  particle  located  at  position  (i,J)  [8]. 

Investigation  into  Accuracy  of  an  Interior  Point 

In  this  example,  the  material  is  elastic  and  does  not  attenuate,  so  the  mass-spring-lattice 

model  (MSLM)  is  used.  The  phase  speed  accuracy  can  be  investigated  using  the  MSLM  and  a 
four-node  FEM  with  an  arbitrarily  oriented  plane  wave  propagating  through  an  unbounded 
elastic  media. 

The  equations  of  motion  in  a  plane  strain  isotropic  medium  having  mass  proportional 
damping  expressed  in  Cartesian  coordinates,  are 


d2u  /  \d°u  /  \  d2v  d2u  p  du 


(4) 


d2v  /  \32v  /  \  d2u  d2v  pdv 

p—  =  (X  +  2P—  +  (X  +  P—  +  IJ—--- 


(5) 


where  u  is  the  displacement  in  the  x-direction,  v  is  the  displacement  in  the  y-di  recti  on,  p  is  the 
density,  r  is  a  time  constant,  and  X  and  p  are  Lame  elastic  constants. 


16 


The  four  node  finite  element  model  discretization  of  eqn.  (4)  and  eqn.  (5)  give  the 
following  equations  written  in  component  fonn  at  a  particle  position  (ij)  at  time  t.  [Appendix  A] 


\  i  \  (  i  i  \  (  i 

t  1  P  1  t—At  1  P  1  1 

U:  ,  p - r-  H - — —  M  •  /? - ; - +2  iq.  .  p - r- 

7 1  (At)  r  2At  J  1,7 1  (At)2  r  2At  J  1,7 1  (At)2 


3  +  p  /  {  ,  t  \  3  ^  It  ,t  ^ I  ) 

+  h2  {  umj+  ui-u  -2  Uu)—^r\  uu+ 1+  uij-i  -2  uuj) 


+  - —  2^  {'ui+hj+ 1  +,Mi-ij+i  +X-i,./-i  +\-+i,./-i  “4'w,  - ) 
n 

4^  +  4^//  t  t  t  \ 

"l  ^2  l  Vi+l,t+l  _  +  V/-l  J-l  _  Vi+1  J-l  j 


(At)  r  2 At 


,+Ai v,  p  — —  +  ^  —  =-'-A/ v,.  p  — —  -  ^  —  +2 ‘v  \p  — 


(At)  t  2 At 


-^T^('v,/+i  +' Aj-i 2v;V )-  ^(‘v/+lj-  + Vij  -2' v, . ) 


2  A  +  2  u  / 1 

+  — 7^ — 1  vi+ij+i+  vmj-\  +  vi- ij-i+  vi-ij+i  ~ 

h 


4'v.J 


\A  +  \ptt  t  t 

+ - ~2 - 1  umj+\~  M/+ij-i+ 

h 


i  Hi-y+i) 


where  At  is  the  numerical  time  step  and  h  is  the  grid  spacing.  (A  similar  analysis  for  the  MSLM 
can  be  found  in  reference  [8].)  For  stability,  the  Courant  number  C  must  satisfy 

Cs£t^<  i  (g 

h 

where  cP  is  the  longitudinal  wave  speed  given  by 

„  _  \A  +  2A 


The  shear  wave  speed  cs  is  given  by 


CS  =J~ 
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Analysis  of  accuracy 

Taking  the  two  dimensional  discrete  Fourier  transform  of  eqn  (3)  and  (4)  and  forming  the 
amplification  equation  yields  [8] 


f+Au=G'u 


(ID 


where 


‘u  =  [*U  'V  A,U  ‘~A,v] 
t+At  u  =  [t+AtU  t+AtV  ‘U  ,v\ 


a 


G  = 


d 

1 

0 


b  c 
e  0 
0  0 
1  0 


0 

/ 

0 

0 


(12) 

(13) 


(14) 


and  where 


vv 


2(2/1 +  3//)n 

3 h2  , 

4(1 +  3//)  A 
6h2  > 


(cos(kv/i)-l) 


'  22 
3  h2 


(cos(kvh'j(cos(kxh))-\) 


P 

(At)2 


+ 


P 

2xAt 


(15) 


-  (sm(kxh)sm(kvh)) 

h_ _ ' 

(At)  It  At 


{(At)2 

It  At  J 

P,+ 

P 

(At)2 

It  At 

(16) 


(17) 
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(18) 


A  +  ju 


d  = 


(s  i  n  (/cv  A)  s  i  n  (A:/; )) 


P  +  P 


(A  t)  2tA  t 


2(2  A  +  3/u) 


\A 


3  h2 
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In  eqns.  (12)  through  (17) 
kx  =  k  cos  6 

k„  =  ksm6 


(20) 


(21) 

(22) 


The  four  eigenvectors  of  the  amplification  matrix  G  are  £+/>,  £-p>  and  Cs,  where  the  subscripts 

+  or  -  refer  to  the  positive  and  negative  phases. 

The  positive  phase  change  in  one  time  step  is  found  from  the  exact  dispersion  relation  as 

(23) 


^exacAt  =  kcpAt 


VexacA*  =  kcSA 


and  the  positive  phase  change  from  the  numerical  approximation  is  [13] 

^numerical**  = 

VnumericaAt  = 


(24) 

(25) 

(26) 


where  fr,  and  its  ore  the  eigenvectors  that  correspond  to  the  pressure  or  P-  and  shear  or  S-waves. 
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which  is  simply  the  percent  error  in  wave  speed. 

Consider  the  exact  and  numerical  dispersion  relation  for  an  elastic  material  with  v=0.30 
and  a  plane  wave  propagating  along  the  x-axis  (0=0)  as  shown  in  Fig.  4.  For  this  example  and  for 
the  remainder  of  this  paper,  C=1  will  be  used.  As  seen  in  the  figure,  the  phase  error  of  P- waves 
of  both  FEM  and  MSLM  rapidly  decrease  despite  the  small  number  of  grid  spaces  per 
wavelength.  This  phenomena  is  due  to  the  fact  that  the  FEM  and  the  MSLM  solutions  are  equal 
to  the  exact  solution  for  this  Poisson’s  ratio  v=0.30  and  angle  of  orientation  0  =0°.  The  S-wave 
propagation  of  the  FEM  and  MSLM  differ  as  the  nonnalized  wavelength  increases.  It  is 
interesting  to  note  that  the  S'-wave  phase  speed  error  is  the  same  for  both  models  at  Poisson’s 
ratio  v  =0.3  and  angle  of  orientation  6  =0°. 

Fig.  5  shows  the  case  where  Poisson’s  ratio  v  =0.3  and  angle  of  orientation  6  =45°. 
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The  only  difference  between  Fig.  4  and  Fig.  5  is  the  angle  of  orientation.  The  FEM  P- wave  error 
is  now  significantly  different  due  to  the  change  in  9.  The  P- wave  solution  changed  from  the 
exact  solution  to  a  solution  that  has  an  error  inversely  proportional  to  the  number  of  grid  spaces. 
In  the  FEM,  the  P-  and  S-  wave  phase  speed  errors  are  nearly  equal  but  the  errors  in  the  P-  and  S- 
waves  for  the  MSLM  are  different  by  a  factor  of  two.  The  unusual  behavior  of  the  MSLM  S- 
wave  is  disregarded  due  to  the  small  number  of  grid  spaces  per  wavelength.  At  two  grid  spaces 
per  wavelength,  both  models  are  extremely  inaccurate. 

While  it  is  well  known  that  the  phase  errors  in  most  numerical  models  vary  as  6  varies,  it 
is  interesting  to  note  that  the  phase  error  in  each  model  varies  as  v  varies  as  well.  In  Fig.  6,9  = 
45°  and  v  =  0.2.  The  phase  error  of  the  MSLM  is  clearly  less  than  the  phase  error  of  the  FEM. 

The  difference  in  error  between  the  FEM  and  the  MSLM  equates  to  roughly  one  less  grid 
spacing  per  wavelength. 

In  Fig.  7,9  =  45°  and  v  =0.4.  The  phase  speed  error  for  the  P- wave  speed  in  the  FEM  is 
slightly  larger  than  the  phase  speed  error  in  the  MSLM.  However,  the  .S'-wavc  speed  error  is 
nearly  one  order  of  magnitude  smaller  in  the  FEM  than  in  the  MSLM. 
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Figure  4.  Percent  error  in  phase  speed  as  a  function  of  grid  spaces  per  wavelength,  where  Poisson’s  ratio  v  =0.3  and 
angle  of  orientation  6  =0°  for  (a)  FEM  and  (b)  MSLM. 
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Figure  5.  Percent  error  in  phase  speed  as  a  function  of  grid  spaces  per  wavelength,  where  the  Poisson’s  ratio  v  =0.3 
and  angle  of  orientation  9  =45°. 


Figure  6.  Percent  error  in  phase  speed  as  a  function  of  grid  spaces  per  wavelength,  where  Poisson’s  ratio  v  =0.2  and 
angle  of  orientation  9  =45°. 
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Figure  7.  Percent  error  in  phase  speed  as  a  function  of  grid  spaces  per  wavelength,  where  Poisson’s  ratio  v  =0.4  and 
angle  of  orientation  9  =45°. 


Typically,  in  engineering  applications  an  acceptable  error  is  known  and  grid  spaces  per 
wavelength  N  are  determined  and  given  by 
2  n 


NP  = 


kph 


AE  = 


2  n 
ksh 


(31) 


(32) 


where  NP  and  N$  are  the  number  of  grid  spaces  per  P  and  S  wavelengths  respectively,  kP  and  ks 
are  the  respective  wavenumbers,  and  h  is  the  grid  space. 

Figure  9  is  the  number  of  grid  spaces  required  by  the  MSLM  to  achieve  1%  error  or  less 
in  phase  speed  as  functions  of  Poisson’s  ratio  and  angle  of  incidence.  Fig.  10  is  the  same  plot  for 
the  FEM.  Note  that  for  both  the  FEM  and  MSLM,  Np  and  Ns  are  symmetric  about  6  =  45°,  which 
follows  from  the  symmetry  of  the  models  for  interior  particles  as  shown  in  Fig.  2  and  Fig.  3.  In 
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both  models,  the  phase  error  of  P- waves  decays  to  zero  as  0  — >0°.  Holding  6  =0°  or  90°,  Ns  is  a 
monotonically  increasing  function  for  both  FEM  and  MSLM 

There  are  differences  in  the  models.  In  the  MSLM,  Np  decreases  as  v  increases  while 
holding  0  =  45°.  In  the  FEM,  holding  0  =45°  and  increasing  v  increases  Np  slightly.  Holding 
#=45°,  Ns  for  MSLM  has  a  minimum  at  v  =0.3  before  increasing  dramatically.  The  FEM  has  a 
minimum  Ns  at  v  =  0.43  before  increasing.  The  errors  in  S-wave  speed  for  the  case  when  0=45° 
and  v  varies  from  0.0  to  0.5  is  shown  in  Fig.  8. 
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Figure  9.  Number  of  grid  spaces  per  wavelength  required  to  achieve  1%  error  or  less  phase  speed  error  in  FEM  as  a  function  of  Poisson’s  ratio  v 
angle  of  incidence  9  with  respect  to  horizontal  axis  for  (a)  P  waves  and  (b)  S  waves. 


Figure  10.  Number  of  grid  spaces  required  for  1%  phase  speed  error  for  S-waves  as  v  varies  from  0.0  to  0.5  angle  of 
incidence  6  =45°. 


Investigation  into  Accuracy  of  a  Point  at  a  Traction  Free  Surface 

Wave  propagation  problems  involving  a  surface  excitation  and  a  surface  response  in  an 

elastic  half-space  have  become  known  as  Lamb’s  problems  due  to  the  efforts  of  Horace  Lamb  in 
1904.  Lamb’s  work  is  based  on  Rayleigh’s  discovery  of  surface  waves  in  1887.  Lamb  discovered 
Rayleigh  waves  are  a  direct  function  of  the  source  kinematics  and  that  P-  and  S-waves  are  a 
function  of  the  time  derivative  of  the  source  function  [14].  An  analytical  solution  exists  for  these 
types  of  problems.  (Refer  to  [8]  or  [15]  for  a  detailed  discussion  of  the  analytical  solution.) 

Previously  in  this  paper,  FEM  and  MSLM  were  compared  to  each  other  for  phase 
speed  accuracy  at  an  interior  point.  FEM  and  MSLM  displacement  accuracy  will  be  compared  to 
an  analytical  solution  for  a  point  on  a  traction-free  boundary. 

Two  dimensional  schematics  of  the  FEM  and  the  MSLM  discretization  of  a  plane  strain 
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elastic  solid  at  a  traction  free  boundary  are  shown  in  Fig.  1 1  and  Fig.  12.  The  corresponding 
equations  of  motion  in  indicial  notation  for  the  FEM  are 


U'-x-2X,+'*X)+£^(- 
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where  At  is  the  numerical  time  step,  and  Uy  and  Vy  are  horizontal  and  vertical  displacements 
respectively.  (A  detailed  analysis  of  the  corresponding  equations  of  motion  of  the  MSLM  can  be 
found  in  reference  [8].) 
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boundary  face 


h 

Figure  11.  Schematic  of  an  FEM  at  surface  plane-strain  particle  located  at  position 

boundary  face 


Figure  12.  Schematic  of  an  MSLM  at  surface  plane-strain  particle  located  at  position  (;',  /)  [8]. 

Input  Signal  Shape 

The  surface  excitations  in  the  two-dimensional  numerical  simulations  are  point  loads. 
The  frequency  content  of  the  time-varying  point  load  is  dictated  by  the  Gaussian-modulated 
cosinusoid  depicted  in  Fig.  13.  The  numerical  function  is 
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u{0,t)  =  up  exp 


(35) 


'  ~\(2  nfj  -  3)2]cos(2^<  -  3  fJJ1 ) 


(Time)/(Period  at  Center  Frequency) 

Figure  13.  Gaussian-modulated  cosinusoidal  input. 

where  u p  is  the  peak  displacement,  fa  is  the  standard  deviation  cyclic  frequency,  and  fc  is  the 

center  cyclic  frequency.  The  maximum  effective  frequency  of  the  input  signal  is  fc  +  3  fa . 

Maximum  effective  frequency  relates  to  the  minimum  propagating  wavelength.  This  minimum 
wavelength  is  used  in  accuracy  calculations. 

Figure  14  shows  the  surface  displacement  of  a  particle  3.5  P- wavelengths  from  the  point 
source  in  the  MSLM.  Figure  15  shows  the  surface  displacement  of  a  particle  3.5  P-wavelengths 
from  the  point  source  in  the  FEM.  Twenty  grid  spaces  per  minimum  wavelength  are  used  in  each 
model.  This  grid  spacing  is  selected  to  ensure  less  than  1%  error  in  phase  speed  for  an  interior 
particle  for  both  the  FEM  and  the  MSLM. 
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Accuracy 

Both  the  FEM  and  MSLM  reproduce  Lamb’s  solution  for  surface  displacements  of  an 
elastic  material  due  to  surface  excitation.  As  shown  in  Fig.  14  and  Fig.  15,  the  difference 
between  the  analytical  solution  and  the  numerical  solution  is  small.  In  the  FEM  and  MSLM,  the 
phase  error  in  both  the  horizontal  and  vertical  directions  is  less  than  2%.  A  detailed  summary  of 
Lamb’s  solution  appears  in  Appendix  C. 


Figure  14.  Horizontal  and  vertical  displacements  of  the  surface  particle  located  3.5  P-wavelengths  from  a  point 
source  in  the  MSLM. 
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Figure  15.  Horizontal  and  vertical  displacements  of  the  surface  particle  located  3.5  P-wavelengths  from  a  point 
source  in  the  FEM  mode. 


Attenuation  analysis  of  FEM  and  MSDLM 

Attenuation  is  a  decrease  in  the  amplitude  of  a  wave  due  to  the  material’s  absorption  of 

energy.  Damping  is  the  tendency  of  a  material  or  system  to  reduce  the  amplitude  of  oscillations. 
The  attenuation  of  a  signal  is  the  result  of  damping  in  a  material  through  which  a  signal  is 
traveling.  In  lumped  parameter  models  damping  is  modeled  by  an  element  retarding  force  is 
proportional  to  the  velocity  across  it.  As  shown  in  Fig.  1,  the  MSDLM  includes  dashpots,  which 
act  as  dampers. 

The  FEM  does  not  model  damping  in  the  above  manner.  A  damping  matrix  cannot  be 
constructed  in  the  same  way  as  the  stiffness  matrix  in  Appendix  A.  Rayleigh  damping  is  used 
instead.  The  equation  for  Rayleigh  damping  is  [12] 

C  =  aM  +  (3K  (36) 
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Where  C  is  the  damping  matrix,  M  is  the  mass  matrix,  K  is  the  stiffness  matrix,  and  a  and  [>  are 
constants  determined  from  two  given  damping  ratios  that  correspond  to  two  unequal  frequencies 
of  vibration.  Thus  damping  in  the  FEM  is  not  proportional  to  velocity  but  is  proportional  to 
mass. 

In  this  thesis,  one  characteristic  of  damping  is  penetration  depth  ( Q )  (Appendix  D).  This 
is  the  depth  an  input  signal  travels  into  a  material  before  the  amplitude  of  the  signal  is  reduced  to 
e (about  4%)  of  its  original  amplitude.  Initial  results  indicate  that  a  damped  FEM  does  agree 
with  numerical  solutions  determined  by  the  MSDLM.  In  Fig.  16,  the  phase  speed  error  of  a 
material  with  Q=  25  is  shown. 
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Figure  16.  Phase  speed  error  of  a  FEM  of  a  material  with  Q= 25  and  angle  of  orientation  8  =45°. 

In  Fig.  17,  the  same  material  is  modeled  with  varying  grid  spaces  per  wavelength.  This 
shows  good  agreement  with  an  analysis  of  the  MSDLM  in  reference  [8]. 
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10 


Figure  17.  Percent  error  in  phase  speed  as  a  function  of  grid  spaces  per  wavelength,  where  Q~L 5,  angle  of 
orientation  0  =45°. 

Computational  cost  of  FEM  and  MSDLM  models 

Accuracy  is  not  the  only  issue  of  importance  for  numerical  simulation.  Computational 

cost  is  another  issue  considered  in  this  study.  Theoretically,  given  enough  time  and  computing 
power,  extremely  accurate  numerical  models  can  be  generated  using  any  numerical  simulation 
method.  In  most  engineering  analyses,  there  are  limits  on  time  or  computing  power  available  for 
a  specific  model.  A  trade-off  study  can  be  completed  in  order  to  maximize  accuracy  while 
minimizing  cost.  In  most  cases,  there  comes  a  point  where  additional  improvement  in  accuracy  is 
not  beneficial  enough  to  justify  the  additional  computational  effort.  This  point  varies  depending 
on  the  type  of  engineering  problem  being  analyzed. 
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One  method  of  measuring  computational  cost  is  by  estimating  the  number  of  floating 
point  operations  (FLOPS)  per  time  step  or  iteration.  For  the  FEM,  the  MSLM  and  the  MSDLM, 
this  involves  detennining  the  number  of  nodes  and  multiplying  that  number  by  the  numerical 
operations  (multiplies,  adds,  divides  and  subtracts)  per  node. 

Since  the  FEM  and  MSLM  use  the  Central  Difference  Method,  the  number  of  FLOPS 
can  be  easily  determined.  For  the  FEM,  the  number  of  numerical  operations  can  be  found  by 
examining  eqn.  (6)  and  eqn.  (7).  The  number  of  nodes  can  be  determined  by  using  the  desired 
accuracy  to  find  the  number  of  elements  needed  to  reach  this  accuracy.  There  are  38  numerical 
operations  per  time  step  per  node  for  the  FEM. 

The  discretized  equations  of  motion  for  the  MSLM  are 
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Figure  18.  Number  of  FLOPS  per  ^-Wavelength  required  for  1%  phase  speed  error  for  S-waves  as  v  varies  from  0.0 
to  0.5  angle  of  incidence  0  =45°. 

Determining  the  FLOPS  for  the  MSDLM  is  more  difficult,  since  the  integration  method 
used  is  a  fourth  order  Runge-Kutta.  Once  the  number  of  numerical  operations  is  determined  by 
examining  eqn.  (36)  and  eqn.  (37),  that  number  must  be  multiplied  by  four.  Some  computational 
savings  are  found  since  the  MSDLM  requires  a  Courant  number  of  1.3.  This  larger  Courant 
number  leads  to  a  larger  time  step,  therefore,  there  are  fewer  nodes  in  the  model.  The  stress- 
dynamic  equations  for  the  MSDLM  are 
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There  are  288  numerical  operations  per  time  step  per  node  for  the  MSDLM.  Stress- 
dynamic  equations  are  used  for  the  MSDLM  to  account  for  changes  in  volumetric  forces.  These 
forces  must  be  included  for  the  stress  relaxation  mechanism  in  the  MSDLM. 
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The  surface  stress-dynamic  equations  for  the  MSDLM  and  the  equations  of  motion  for 
the  FEM  and  the  MSLM  are  different  from  the  interior  equations.  For  most  engineering 
problems,  the  number  of  surface  particles  is  minuscule  compared  to  the  number  of  interior 
particles  and  so  are  neglected  in  this  study.  The  computational  cost  of  a  four-node  FEM  is 
approximately  the  same  as  that  for  the  MSLM.  The  computational  cost  of  the  MSDLM  is  on  the 
order  of  5  times  greater  than  the  cost  of  the  FEM  and  MSLM. 

Conclusions 

Finite  element,  mass-spring  lattice,  and  mass-spring-dashpot  lattice  models  are  powerful 
numerical  simulation  tools.  For  modeling  elastic  materials  with  Poisson’s  ratio  between  0.0  and 
0.2,  the  mass-spring  lattice  model  has  the  lowest  computational  cost  for  phase  speed  error  less 
than  1%.  For  modeling  elastic  materials  with  Poisson’s  ratio  between  0.35  and  0.45,  the  finite 
element  model  is  more  accurate  but  costs  just  slightly  more  than  the  mass-spring  lattice  model. 
For  modeling  materials  with  attenuative  properties,  the  MSDLM  is  more  accurate  but  is  nearly  5 
times  more  expensive  than  the  finite  element  model.  Both  the  FEM  and  MSDLM  model 
attenuative  materials  accurately. 
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APPENDIX  A  -  Derivation  of  the  indicial  notation  for  the  4-node  finite 
element  model 

The  finite  element  method  (FEM)  is  a  widely  used  tool  in  engineering  analysis  [A.  1] .  The 
four-node  element  is  one  of  the  simplest,  but  it  is  still  a  powerful  and  useful  tool  in  the  solution 
of  practical  engineering  problems.  In  this  appendix,  a  plane  strain  isotropic  medium  having 
mass-  proportional  damping  is  discretized. 

The  equations  of  motion  in  the  interior  of  a  plane  strain  isotropic  medium  having  mass- 
proportional  damping  is  expressed  in  Cartesian  coordinates  as 


d°u  \d2u  / .  \  d2v  d2u  p  du 


(A.l) 


82v  /,  n.S2v  /,  \  d 2 u  82v  p  dv 


( A.2 ) 


where  p  is  density,  u  and  v  are  horizontal  and  vertical  displacements,  respectively,  X  and  p  are 
Lame  elastic  constants  and  r  is  a  characteristic  time. 

The  corresponding  elemental  equilibrium  equation  for  a  lumped  mass  element  can  be 
expressed  in  matrix  form  as 


f+A/  (e/)  V  (e/)  ,  f-Af  (e/)  i 

M(e/)  u  z  u  +  u  |  i  M(g/|  u _ U  I  K(e/)?u(e/)  =  F(e7> 

'2  -  2A  t 


t+M  (el)  t-At  (el) 


(a  ty 


where  M'cl>  is  the  elemental  mass  matrix. 


(A. 3) 


(el)  _  Ph 


1VT  = 


(A.4) 
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and  'uid)  is  the  elemental  displacement  vector 


A  schematic  of  a  square  four-node  element  having  grid  spacing  h  is  shown  in  Figure  Al. 


Figure  A  1.  Schematic  of  four-node  element. 


The  elemental  stiffness  matrix  K'Wy  may  be  calculated  as 


h /  y 

K,el)  =  J  | B <el) T C <eI) B (e,) dxdy  (A.6) 

-h/  -h/ 

72  72 
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n(el)  is  the  displacement  interpolation  matrix  for  a  four-node  square  element  i  with  sides  of 
length  h 


H(el) 


hx  h2  h3  h4  0  0  0  0 

0  0  0  0  A,  h2  h2  h4 


(A.  7) 


where 


fh  = 
h2  = 

K  = 
h4  = 


i3 


v 


i+ 


A 


2x 


V 


1 - 

h  A 


1  + 


1-^ 


v 


l- 


A 


1  +  — 


V 


1- 


A 


2y^ 

h  y 

Vi 
h  , 

V 

h  , 

V, 

2hy 


(A. 8) 

(A. 9) 

(A.  10) 

(A.  11) 


B (el)  is  the  corresponding  strain-displacement  matrix  corresponding  to  the  local  element  degrees 
of  freedom,  given  as 


-sf 

i _ 

Kx 

K 

Kx 

0 

0 

0 

0 

B(el)  = 

0 

0 

0 

0 

Ky 

Ky 

Ky 

Ky 

Ky 

Ky 

Ky 

Ky 

Kx 

Kx 

Kx 

Kx 

where 


(A.  12) 
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C (eI>  is  the  generalized  stress-strain  matrix  for  iso 


conditions,  given  as 


C(e,) 


A  +  2/j  A  0 

A  A  +  2ju  0 

0  0  ju 


(A.  13) 


(A.  14) 


(A.  15) 


(A.  16) 


(A.  17) 


(A.  18) 


(A.  19) 


(A. 21) 


The  elemental  K  matrix  is 


K(e,) 


K  k2  k3  k4  k5  k6  k7  k& 

kx  k4  k3  k&  k7  k6  k5 

kx  k2  k7  ks  k5  k6 

K  k6  k5  ks  k7 

kx  k4  k3  k2 

k  |  k2  k3 

kx  k4 
kx 


where 


k\  —  —  A,  +  u 
3 


,  "I;  1 

K2  =  - A - U 

3  2 


Z  — 1  ;  1 

h=  —A--fU 

6  2 


k4=  —A 
6 


,  1,1 

=  —A  +  —  u 
4  4 


,  _  1  ;  1 

—  —  A - u 

4  4 


z  -  1  ;  1 

K'] - A - LI 

4  4 


(A.22) 


(A. 23) 

(A. 24) 

(A. 25) 

(A. 26) 

(A. 27) 

(A. 28) 

(A. 29) 
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(A. 30) 


k% - A  — u 

4  4  ^ 

Interior  Element 

Through  careful  examination  of  the  matrices  centered  at  (ij),  the  global  degrees  of 
freedom  can  be  identified  through  use  of  the  local  degrees  of  freedom  and  the  matrices  listed 
above. 


i-l,j-l  i,j+l  i+l,j+l 


Figure  A  2.  Interior  elemental  matrices  centered  at  position  (ij). 


48 


The  four-node  finite  element  discretization  of  eqns.  (A.l)  and  (A.  2)  yields  the  following 
equations  written  in  component  fonn  at  position  (/,/)  and  time  t. 


-2'“U+'*X)=^(WV,J  ~2'UJ 


('Uij+i+'uij-i  -2\J 
A  +  jij  ( 


f  A  +  \  LI  j 

+  - - ~~ — I'm..,,  u;  ,  u:  ,  ,+' U; 

h 


1  (+1,7+1  1  i-1,7+1  1  “i-lj-l  1  "i+lj-l  ^  “ij 


4'"J 


4"  1  4<“  r..  ,  t 


2  \  Vi+\,j+\  Vi-\,j+ 1  Vi+l,j-l , 


h 


P  1  f  1-A1  ,  /+A1  ) 

- 1-  M,  ,+  M,  ,1 

r  2AA 


(A. 31) 


A  +  fi 


—  ('“AV  .  -2V  .  +  '+AV  .)  = 

At21  A2 


v  “v/./  v'yj-  1I^1('vi+'vi-'2vu) 


,2 


+ 


A  v 
+  \/u 


+ 


h  v 

+  \n  (t 

h2  ' 

P  1  /  /-A 

r  2At V 


V-  ,  . 

-2'v, 

>-1.2 

hj  ) 

’i+U+i 

+  'Vi  +  1.2-l 

t 

( 1  +  1, 2  +  1 

~  U  1  +  1.2“ 1 

Vi,2  + 

f+Af  ) 

ViJ 

■'VHM  +  'Wi  ~4'Vu) 


(A. 32) 


where  At  is  the  numerical  time  step. 
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Surface  Element 


Through  examination  of  the  matrices  centered  at  (ij),  the  global  degrees  of  freedom  can 
be  identified  through  use  of  the  local  degrees  of  freedom.  A  two-dimensional  schematic  of  the 
FEM  discretization  of  a  plane  strain  elastic  solid  along  a  traction  boundary  is  shown  in  Figure  A 
3. 


boundary  face 


Figure  A  3.  Schematic  of  a  surface  element  centered  at  (i,j). 


The  corresponding  equations  of  motion  are 
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j^ft  t  ) 

\X  +  H(t  , t 

+ — 77 — l  “i-u-l+  “<+i.y 


jA+pft 


z  2A t 


2'“«) 


+  li±lii('v  _'v 

+  ^2  V  VMJ-1  Vi+l,y-l 


J-l  2'"/./) 


)-^  W'y.  -'v.  ' 

/  7,2  V  '-1,2  '+0, 


(A. 33) 


—  ('-A'v.  .  -2'v.  ,+'+V  ,)+^  — (-'-V  .+,+A,v.  ■)=  3  ^  +  (<  '  .) 

V  ',2  ',2  ',2/  ,  1.3  '.2  ',2/  ,2  V  ',2-1  '-2/ 


r  2At 


y/l  +  2//  / 
/22  ^ 


^('Vl-U+'Vwj  -2^', 2) 

7A  +  V(t  ,t  V  ) 

-Tl—l  Vu-1+  vi+u_i  -2  V,.j 


,,2  l  W<-1,2-1  M '+1,2-1/  ^2  l  w'+l,2 


M'+l,2  Ui~\,j  t 


(A. 34) 


Where  A t  is  the  numerical  time  step,  and  u,i7  and  Vy  are  the  horizontal  and  vertical  displacements, 
respectively. 
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APPENDIX  B  -  Dispersion  Relation/  High  frequency  assumption 

In  order  to  make  numerical  analysis  of  elastic  wave  phenomena  possible,  assumptions 

concerning  the  dispersion  relation  and  frequency  of  the  wave  must  be  made.  In  this  appendix,  the 
dispersion  relation  for  a  material  having  mass-proportional  damping  is  derived. 

Equilibrium  equation 

The  equation  of  motion  for  a  one-dimensional  material  having  mass-proportional 
damping  is 


pd2u/  i _P_du 
7dx2+  r 


It/ 

2  +  /dx~p  /dt2 


(B.l) 


where  E  is  the  governing  elastic  constant,  u  is  the  displacement,  p  is  the  density,  and  r  is  a 
characteristic  time. 

Consider  a  propagating  harmonic  wave  having  the  form 
u(x,t)  = 


(B-2) 


where  U0  is  the  peak  displacement,  a  is  the  attenuation,  k  is  the  wavenumber,  and  co  is  the  radial 
frequency. 

Solving  for  the  partial  derivatives  yields 


32m/_2  =U0(a2 -k2)e-ax+m-,ot)  - U u (2aik)e~ax+Ukx~0X) 

(B.3) 

du/^  =  u  oict)e~ax+ilkx~,u) 

(B-4) 

d2u /  -  _7T  p-tvc+'Va-M) 

Vdt2~  UoCOe 

(B.5) 

After  substituting  eqns.  (B.3)  through  (B.5)  into  eqn.  (B.l)  and  simplifying,  the 
dispersion  relation  becomes 
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E(a?  -  k2)  +  pco1  =  E(2aik)  -  —  ico 

X 

Separating  the  real  and  imaginary  parts  of  eqn.  (B.6)  yields 

(B.6) 

E(a2  -k2)  +  pco2  =  0 

(B.7) 

E{2aik)-—ico  =  0 
r 

Solving  eqn.  (B.8)  for  the  wavenumber  yields 

(B.8) 

k  -  _p^_ 
xE2a 


(B-9) 


Substituting  eqn.  (B.9)  into  eqn.  (B.7)  yields 


a  - 


r  _ 


pco 

xE2a 


+  co  —  0 


p 


(B.10) 


Let 


~  =  C2  max  (B.ll) 

P 

Multipyling  eqn  (B.10)  by  a  and  substituting  eqn.  (B.l  1)  into  eqn.  (B.10)  yields 


2  2  2 

4  co  co  a 
a  - -  + - ^  =  0 


r  4 c„ 


4 


(B.12) 
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Solving  for  a2  yields 


(B-13) 

(B.14) 

(B.15) 

(B.16) 

(B.17) 
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APPENDIX  C  -  Lamb’s  Solution 

This  appendix  summarizes  Lamb’s  solution  for  surface  displacements  due  to  surface 
excitation  [C- 1  ] .  Consider  a  normal  point  line  source,  having  peak  magnitude  Q  and  temporal 
variation  q(t),  acting  at  the  origin  of  a  half-space  (v<0)  having  Lame  constants  X  and  p,  and 
density  p. 


For  x  <  0 ,  the  horizontal  and  vertical  surface  displacements,  u  and  v  respectively,  are 


— Hq(t  -  sRx) - J 

//  II  7T  " 


Qult  \  02}  sr2e{id2  -sT2\le2  -sL2  ^sT2  -e2 


M  n  {id2  -  ST2  )4  + 1 694  [o2  -  s2  )(5r2  -  92 ) 


-  (k)dd  (C.l) 


where  sT  is  the  transverse  wave  slowness  (inverse  wave  velocity)  given  by 


sT  = 

s,  is  the  longitudinal  wave  slowness  given  by 


s 


sR  is  the  Rayleigh  wave  slowness  given  by  the  real  root  of  the  equation 


F 


(C.3) 


(C.4) 


(C.5) 


(C.6) 
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where 


4|  -  2s j1  sT2 s R  +3bi2  +sT2 


<N 

2 

/  2 

2  1 

JSR 

~SL  ' 

\JSR 

—  sT 

y 

1  2" 

2 

(N 

<N 

UJ  rp  ij  p 

V SR 

~SL  ' 

Vs  R 

—  sT 

F  = 


P  in  eqn.  (C.2)  denotes  the  principal  value  of  the  integral  [C-2].  A  non-integrable  singularity 
exists  at  9=  sR  and  the  integral  must  be  defined  as 


p | /fete = iim  j  /fete  + }  /fete 


(C.8) 


V  a 


where  a<c<b  and  the  non-integrable  singularity  exists  at  f(c). 

The  surface  displacements  for  x<0  are  given  by  replacing  (/  -  fx)  with  (/  +  fx)  in  eqns. 
(C.l)  and  (C.2)  and  reversing  the  sign  of  the  horizontal  displacements  in  eqn.  (C.l). 
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APPENDIX  D  -  Penetration  Depth 

Let  the  penetration  depth  Q  of  an  attenuative  material  be  the  number  of  wavelengths  required  for 
a  cyclic  input  signal  to  decrease  in  amplitude  to  e where  a  is  the  attenuation  and  x  is  the 
distance. 


e~aQl  =  e~* 

(D.l) 

where  X  is  the  wavelength 

k; 

ll 

Oj 

(D.2) 

Q  =  — 

aX 

(D.3) 

For  a  standard  linear  solid  [Dl]  having  a  relaxation  time  z  ,  under  the  frequency  assumption 

on  »1  (D.4) 

where  co  is  the  circular  frequency.  Frequency-independent  attenuation  is 

a  =  —  (D.5) 

2  TC 

where  c  is  maximum  phase  speed.  Recall, 

c  =  W  (D.6) 

where/is  cyclic  frequency.  The  relationship  between  co  and/is 

co  =  2nf  (D.7) 

therefore 


/  = 


co 

In 


Substituting  eqn.  (D-5)  into  eqn.  (D-3)  yields 
nlzc 


1 


=  Q 


Substituting  eqn.  (D-6)  into  eqn.  (D-9)  yields 


(D.8) 


(D.9) 
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nlrf  =  Q 


(D.10) 


Substituting  eqn.  (D-8)  into  (D-10)  yields 
to)  =  Q 

Recall, 

co 

c  =  — 
k 

where  k  is  the  wavenumber,  k  can  be  rewritten  as 


(D.l  1) 


(D.12) 


(D.13) 


The  phase  speed  error  between  the  exact  solution  and  the  numerical  approximation  can  be 
written  as 


C  ~  C  max 

error  = - - 

r 

max 


Substituting  eqn.  (D-12)  into  (D-14)  yields 


co 

- c 

i  max 

k 

error  - - 

^  max 


Substituting  eqn.  (D-13)  into  (D-15)  yields 


f  I - T” 

error  =  V2]  1  +  1 1  H - - — - 

[  V  co  T 


(D.14) 


(D.15) 


(D.16) 


Figure  D  1  is  a  plot  of  the  error  with  respect  to  cot  .  The  error  decreases  as  the  frequency 
increases.  As  cot  increases,  Q  increases.  Therefore  a  material  having  a  high  penetration  depth, 
which  is  equivalent  to  a  low  attenuation,  is  modeled  more  accurately. 
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Figure  D  1.  Phase  speed  error  as  a  function  of  cor. 
Error  in  attenuation 
Attenuation  may  be  written  as 


a  = 


co 


yfll 


\/2 


—  1  +  Jl  + 


max  \ 


2  2 

CO  T 


recall  that  numerically 


2  zc 


The  error  between  the  exact  and  numerical  attenuation  can  be  written  as 


(D.17) 


(D.18) 
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error  = 


«max  -  a 


a„ 


(D.19) 


Substituting  eqn.  D-18  and  C-19  into  eqn.  D-17  yields 


error 


=  1-V2z 


\lA 


—  1  +  „  1  + 


2  2 

CO  z 


(D.20) 


Figure  D  2  is  a  plot  of  the  attenuation  error  as  a  function  of  cor  .  As  cot  increases,  the  error 
decreases. 


Figure  D  2.  Attenuation  error  as  a  function  of  an. 
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If  a  material  is  attenuative,  its  penetration  depth  decreases.  Since  penetration  depth  is  a 
function  of  coz ,  the  frequency  decreases  as  well.  This  decrease  in  frequency  causes  errors  in 
phase  speed  modeling. 
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